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1,0  INTRODUCTION 


One  of  the  common  problems  in  numerical  simulations  of  nonlinear 
response  in  an  infinite  domain  is  the  artificial  reflections  that  are 
generated  at  the  finite  boundaries  of  the  computational  grid.  These 
artificial  reflections  are  undesirable  because  they  propagate  throughout 
the  computational  grid,  causing  errors  in  the  computed  response.  The 
arrival  of  the  artificial  reflections  can  be  delayed  by  placing  the  grid 
boundaries  in  a  quiescent  region  that  is  far  from  the  region  of  interest. 
However,  this  approach  is  not  economical  because  the  computational  region 
becomes  quite  large. 


The  general  objective  of  this  effort  is  to  develop  a  transmitting  or 
energy  absorbing  boundary  for  calculations  of  transient,  nonlinear 
response  in  an  infinite  domain.  A  transmitting  or  energy  absorbing 
boundary  attempts  to  eliminate  the  artificial  wave  reflections  from  the 
finite  boundaries  of  the  computational  grid,  allowing  the  computational 
domain  to  be  centered  on  the  region  of  interest. 


A  number  of  nonreflecting  boundaries  have  been  proposed  in  the 
literature.  One  particular  method,  called  the  Incremental  Superposition 
Boundary  (ISB),  appears  to  be  very  efficient  and  computationally  robust 
for  transient  calculations.  The  ISB  method  (Cundall,  et  al.)  is  based  on 
superposition  of  the  solutions  for  a  fixed  boundary  and  for  a  free 
boundary.  Simple  acoustic  theory  shows  that  the  reflected  waves  from  a 
fixed  and  a  free  boundary  are  equal  in  amplitude  but  opposite  in  sign; 
superposition  of  these  two  solutions  will  cancel  the  reflected  waves, 
leaving  the  incident  pulse. 


The  specific  objective  of  the  present  study  is  to  evaluate  the  ISB 
method  as  a  transmitting  or  energy  absorbing  boundary  for  transient  cal- 
culations.  Preliminary  one-dimensional  computational  tests  showed  that 


culations.  Preliminary  one-dimensional  computational  tests  showed  that 
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the  ISB  method  did  not  eliminate  a  reflected  wave  completely;  the  ampli-  Jiced 
tude  of  the  reflected  wave  after  cancellation  was  on  the  order  of  5-12%  of 
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the  amplitude  of  the  incident  wave.  Detailed  analyses  of  the  reflection 
and  cancellation  process  were  performed  to  determine  the  cause  of  the 
errors. 

One  potential  cause  of  the  error  in  the  cancellation  process  is  the 
finite  displacement  of  points  in  the  computational  algorithms.  That  is, 
the  different  displacements  of  a  compression  wave  and  an  expansion  wave 
will  shift  the  pulses,  causing  an  error  in  the  cancellation  process.  A 
simple  perturbation  analysis  shows  that,  for  small  strain  elastic 
response,  the  effect  of  finite  displacement  is  negligible  compared  to  the 
observed  cancellation  error. 

An  analytic  model  was  then  used  to  identify  the  cause  of  the  error  in 
the  cancellation  process  of  the  IS6  method.  This  model,  which  is  based  on 
a  unit  step  pulse  propagating  through  a  one-dimensional  computational 
grid,  shows  that  the  cancellation  is  perfect  when  the  pulse  first  arrives 
at  the  fixed  and  the  free  boundaries;  but  at  the  next  and  any  subsequent 
computational  cycles,  the  cancellation  has  a  substantial  error  if  the 
velocity  at  the  free  boundary  is  not  reset  properly. 

The  analytic  model  also  shows  that  the  appropriate  velocity  at  the 
free  boundary  has  to  be  modified  when  the  computational  time  step  is  less 
than  the  maximum  stable  time  step.  First-  and  second-order  corrections  to 
the  appropriate  velocity  boundary  conditions  have  been  developed  for  time 
steps  less  than  the  maximum  stable  time  step.  Extension  of  the  first-  and 
second-order  corrections  to  two-dimensional  problems  has  not  yet  been 
developed. 

The  theory  and  the  computer  implementation  for  the  basic  ISB  method 
is  discussed  in  Sections  2.1  and  2.2,  respectively.  The  ISB  was  incor¬ 
porated  in  the  Lagrangian,  explicit-in-time,  finite-difference  code 
STEALTH  (Hofmann,  1981),  and  tested  by  performing  calculations  with  a 
cosine  pressure  pulse  propagating  through  a  one-dimensional  grid.  These 
calculations  will  be  presented  in  Section  2.3.  The  effect  of  the  finite 
displacement  of  grid  points  on  the  cancellation  errors  is  discussed  in 
Section  2.4. 


Section  3  presents  the  analytic  model  for  identifying  the  cause  of 
the  cancellation  error  and  the  derivation  of  the  first-  and  second-order 
corrections  for  the  ISB  method.  The  first-  and  second-order  corrections 
for  the  modified  ISB  method  are  then  verified  by  performing  one-dimen¬ 
sional  test  calculations  with  STEALTH.  These  test  calculations  are  also 
presented  in  Section  3. 

Section  4  is  the  summary  and  conclusions  from  this  research. 
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2.0  DESCRIPTION  OF  THE  ISB  METHOD 
2.1  THEORY  OF  THE  ISB  METHOD 

The  concept  for  the  Incremental  Superposition  Boundary  (ISB)  method 
is  based  on  a  non -reflecting  boundary  that  was  developed  by  Smith  (1974) 
for  wave  propagation  problems.  Smith's  method  calculates  the  dynamic 
response  of  a  system  in  an  infinite  domain  by  superimposing  the  complete 
solutions  for  a  fixed  (Dirichlet)  boundary  problem  and  a  free  (Neumann) 
boundary  problem.  The  reflected  waves  from  a  fixed  and  a  free  boundary 
are  equal  in  amplitude  but  opposite  in  sign,  so  superposition  will  cancel 
the  reflected  waves,  leaving  the  incident  wave. 

This  cancellation  can  be  illustrated  for  one-dimensional  acoustic 
waves.  The  displacement,  U,  of  the  incident  wave  in  an  infinite  domain 


may  be  expressed  as: 
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where  A  is  the  amplitude  of  the  wave,  u  is  the  frequency,  c  is  the 
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The  average  of  Equations  (2-4)  and  (2-5)  is  then  the  incident  wave,  Equa¬ 
tion  (2-1).  Similar  arguments  also  apply  to  stresses  and  velocities  (in  a 
linear  medium).  Therefore,  the  solution  for  an  infinite  domain  can  be 
obtained  by  superimposing  and  averaging  the  solutions  for  a  fixed  and  a 
free  boundary  problem. 

The  disadvantage  of  Smith's  method  is  economics:  a  number  of  complete 
solutions  must  be  computed  for  various  combinations  of  boundary  conditions 
and  wave  reflections,  and  is  frequently  simpler  to  just  enlarge  the  com¬ 
putational  domain.  The  ISB  method  avoids  the  calculation  of  complete 
multiple  solutions  through  an  incremental  approach.  This  incremental 
approach  superimposes  the  fixed  and  free  boundary  solutions  in  two  small 
buffer  regions  near  the  boundary,  eliminating  the  reflections  from  the 
boundaries  as  they  occur.  The  dual  calculations  are  only  required  for  the 
small  buffer  regions  with  the  ISB  method. 

2.2  COMPUTER  IMPLEMENTATION  OF  THE  ISB  METHOD 

Figure  2-1  presents  a  schematic  diagram  for  the  structure  of  the 
buffer  regions  and  the  main  grid  for  the  ISB  method.  Two  overlapping 
buffer  regions,  A  and  B,  each  of  three  or  four  zones,  are  connected  inde-  ^ 
pendently  to  the  main  calculation  grid.  Region  A  has  a  fixed  boundary, 
and  region  B  has  a  free  boundary.  For  simplicity,  a  one-dimensional  grid 
is  presented  in  the  diagram;  the  ISB  method  can  be  applied  to  two-dimen¬ 
sional  problems. 

A  wave  that  propagates  from  the  main  grid  enters  the  two  buffer 
regions  simultaneously  at  point  p.  This  wave  will  then  reflect  off  the 
fixed  and  free  boundaries  of  regions  A  and  B.  These  reflected  waves  are 
eliminated  by  replacing  the  existing  variables  in  regions  A  and  B  with  the 
average  values  of  variables  in  the  two  buffer  regions.  (Each  variable  of  a 
zone  in  region  A  and  of  the  corresponding  zone  in  region  B  are  summed  and 
divided  by  two.  The  variables  at  the  nodes  connected  to  the  main  grid  are 
not  averaged.) 
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This  averaging  procedure  does  not  have  to  be  performed  every  cycle 
because  of  the  stability  criterion  for  explicit-in-time  calculations.  The 
ISB  method  is  designed  for  transient  nonlinear  problems.  These  nonlinear 
problems  are  usually  computed  with  an  explicit-in-time  numerical  integra¬ 
tion  scheme.  The  stability  criterion  for  an  explicit  calculation  requires 
that  the  time  step  be  smaller  than  the  time  for  a  wave  to  travel  one  zone 
width.  The  averaging  procedure  is  then  performed  once  every  three  or  four 
cycles  because  the  reflected  waves  cannot  propagate  into  the  main  grid  (or 
to  point  p)  in  less  than  five  cycles. 


fixed 


Direction  of  Wave  Propagation  / 


Buffer  Region  A 


Main  Grid 


Buffer  Region  B 


Figure  2-1  A  Schematic  Diagram  of  the  Buffer  Region  and  the  Main  Grid 
for  the  Incremental  Superposition  Boundary. 


2.3  ONE-DIMENSIONAL  COMPUTATIONAL  TEST  PROBLEM  FOR  THE  ISB  METHOD 


The  ISB  model  has  been  Incorporated  in  the  Lagrangian,  expllcit-in- 
time,  finite-difference  code,  STEALTH  (Hofmann,  1984).  One-dimensional 
calculations  with  a  cosine  pressure  pulse  that  propagates  through  the  grid 
are  performed  to  verify  the  ISB  method.  The  lengths  of  the  main  grid  and 
the  buffer  regions  for  these  test  calculations  are  42  cm  and  4  cm,  respec¬ 
tively.  All  zones  are  one  cm  in  length,  so  the  main  grid  contains  42 
zones  and  each  buffer  region  contains  4  zones. 

Linear  elastic  material  properties  are  assumed  for  the  test  calculat¬ 
ions.  The  equation  of  state  can  be  expressed  as: 

P  -  Po  *  K  -  1) 

where  P  -  pressure  (Mbar), 

Pq  -  initial  pressure  (MBar). 

K  -  bulk  modulus  (MBar), 

Pq  -  reference  density  (gm/cc),  and 
p  -  density  (gm/cc). 

The  initial  pressure,  the  bulk  modulus,  and  the  reference  density  are 
zero,  0.0022  Mbar,  and  1.0  gm/cc,  respectively.  The  pressure  boundary 
condition  that  generates  the  cosine  pressure  pulse  is  given  by: 

P  -  5x10'®  (1  -  cos  ut), 

where  t  -  time  in  fis, 
and  u  -  0.02  t  ns~^. 

Three  calculations  were  performed  with  time  step  safety  factors  of 
0.95,  0.67  and  0.50.  The  time  step  safety  factor  is  the  ratio  of  the 
computational  time  step  to  the  maximum  stable  time  step.  Typical  values 
of  the  time  step  safety  factor  range  from  0.5  to  0.95;  the  default  value 
for  STEALTH  is  0.67. 


Figures  2-2(a),  2'3(a),  and  2-4{a)  present  pressure  time-history 
plots  at  a  point  in  the  main  grid,  3  cm  from  the  inside  edge  of  the  buffer 
region.  The  large  pulse,  centered  at  300  Msec,  is  the  incident  wave;  the 
smaller  waves  after  360  Msec  are  the  reflected  waves  from  the  ISB.  (If 
the  ISB  were  perfect,  only  the  large  pulse  would  appear  in  these  figures.) 
The  three  figures  are  for  time  step  safety  factors  of  0.95,  0.67,  and 
0.50,  respectively. 

Figures  2-2(b),  2'3(b),  and  2'4(b)  present  pressure-distance  plots  at 
500  Msec,  when  the  incident  wave  has  propagated  beyond  the  ISB,  and  the 
rv’flected  wave  is  centered  in  the  main  grid.  (If  the  ISB  were  perfect,  no 
reflected  pulse  would  appear  in  the  figures.)  These  plots  show  the  ampli¬ 
tude  and  time  dependence  of  the  reflected  waves  on  the  same  scale  as  the 
first  set  of  figures.  Figures  2-5  through  2-7  present  the  same  pressure- 
distance  plots  in  an  enlarged  scale. 

These  results  show  that  the  original  ISB  method  does  not  eliminate 
the  reflected  wave  completely;  the  amplitude  of  the  reflected  wave  after 
cancellation  is  on  the  o«^der  of  5-12%  of  the  amplitude  of  the  incident 
wave,  and  is  a  function  of  the  magnitude  of  the  time  step  (or  the  time 
step  safety  factor). 

2.4  EFFECT  OF  FINITE  DISPLACEMENT  OF  THE  COMPUTATIONAL  GRID  POINTS 

Computational  results  from  the  previous  section  show  that  the  ISB 
method  has  a  substantial  error  in  the  cancellation  process.  One  potent¬ 
ial  cause  of  this  cancellation  error  is  the  finite  displacement  of  grid 
points  in  the  computational  algorithms.  That  is,  the  different  displace¬ 
ments  of  a  compression  wave  and  an  expansion  wave  will  shift  the  pulses, 
causing  an  error  in  the  cancellation  process. 
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Figure  2-2  Pressure 


Time  Step  Safety  Factor  is  0.67. 
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Figure  2-5  Pressure-distance  Plot  at  500  ms  for  the  Reflected  Wave 
for  the  Original  ISB.  Time  Step  Safety  Factor  is  0.95. 
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6  Pressure-distance  Plot  at  500  ms  for  the  Reflected  Wave 
for  the  Original  ISB.  Time  Step  Safety  Factor  is  0.67 
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Pressure-distance  Plot  at  500  ms  for  the  Reflected  Wave 
for  the  Original  ISB.  Time  Step  Safety  Factor  is  0.50 


A  simple  perturbation  analysis  can  estimate  the  magnitude  of  the  error 
due  to  the  finite  displacement  of  the  computational  grid.  The  velocities 
of  the  reflected  wave  from  the  fixed  and  the  free  boundary  can  be 
expressed  as: 


*^fixed 


^  (-X  -  ct) 


^  {-X  -  ct) 


(2-6) 


(2-7) 


where  V^.  -  velocity  of  the  reflected  wave  {cm//is), 

A  -  amplitutde  of  the  incident  (or  the  reflected)  wave 
(cm/MS), 

X  *  position  (cm) 

u  -  frequency  (MHz), 

c  -  sound  speed  (cm//is),  and 

t  »  time  (/is) . 

The  error,  E,  of  the  velocity  of  the  reflected  wave  due  to  the  finite 
displacement  of  the  grid  points  can  be  obtained  by  perturbing  Equation 
(2-7)  with  respect  to  u  and  x.  Then  the  error  can  be  expressed  as: 


dx  + 


(2-8) 


Substituting  Equation  (2-7)  into  Equation  (2-8),  the  error,  E,  becomes 


iA  rc 
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(wdx  +  (x  +  ct)dw 


^  1  A  (wdx  +  (x  +  ct)dw)  /  c  \.  (2-9) 

The  maximum  displacement,  dx,  of  a  grid  point  for  the  test  problem  dis¬ 
cussed  in  Section  2-3  is  found  by  integration  of  the  velocity  for  the 
period  of  the  incident  wave.  Acoustic  theory  (Currie,  1974)  shows  that 
the  velocity,  v(t),  of  an  acoustic  wave  is  given  by  : 


v(t)  -  Pb(t)  c  /  K, 


(2-10) 
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where  P|3(t)  is  the  pressure  at  the  boundary,  c  is  the  sound  speed,  and  K 
is  the  bulk  modulus.  The  maximum  displacement  is  then: 


dx 


1 

0 


c  5x10’^ 
K 


(1  -  cos  ut)  dt 


-  0.00337  cm. 


where  T  -  period  of  the  wave  -  100  /is, 
c  -  J{K/p)  -  0.1483  cm/iis, 
p  »  density  »  1  gm/cc,  and 
K  -  bulk  modulus  -  0.022  Mbar. 


The  perturbation  in  u,  du,  can  be  estimated  as: 
.. _ 2ir 


do) 


T  +  dx/c 


=  1. 43x10"^ 


where  u  -  frequency  «  0.02ir  MHz. 

Substituting  these  values  of  dx  and  dw  into  Equation  (2-9),  the  normalized 
error,  1e/a|,  for  the  test  problem  in  Section  2-3  is  found  to  be  0.00286. 
The  same  argument  can  be  applied  for  pressure  and  displacement. 


s 


5 
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The  result  from  this  perturbation  analysis  shows  that  the  cancel¬ 
lation  error  from  the  finite  displacement  is  approximately  0.3%,  well 
below  the  observed  error  and  below  the  level  of  numerical  noise  in  the 
calculation.  (The  value  of  0.3%  is  valid  for  the  specific  pulse  and  mate¬ 
rial  properties  that  are  used  in  the  numerical  calculations.  The  mag¬ 
nitude  of  the  error  will  vary  with  pulse  amplitude,  pulse  shape,  and  bulk 
modulus.) 
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3.0  ANALYSIS  OF  THE  ISB  METHOD 


3.1  IDENTIFICATION  OF  THE  CAUSE  OF  THE  CANCELLATION  ERROR 


The  cause  of  the  cancellation  error  in  the  ISB  method  has  been 
identified  by  a  simple  analytic  model  with  a  step  pulse  propagating 
through  a  computational  grid.  The  pressure  pulse  is  allowed  to  propagate 
into  two  separate  buffer  grids  with  a  fixed  and  a  free  boundary,  and  the 
solutions  of  these  two  grids  are  averaged  at  the  end  of  each  computational 
cycle.  The  initial  pressure  is  po»  and  the  step  pulse  is  generated  by  a 
constant  pressure,  p^,  at  the  upstream  boundary.  This  model  assumes  one¬ 
dimensional  grids  with  two  equal  zones  of  length  L  and  a  linear  elastic 
medium. 


The  equation  of  state  for  a  linear  elastic  material  can  be  expressed 


(3-1) 


where  p  ■  pressure, 

Po  -  initial  pressure, 

K  -  bulk  modulus, 

Pq  ■  reference  density  and, 
p  -  density. 


For  a  small  strain  problem,  the  equations  of  state  and  motion  at  cycle  n 
can  be  expressed  in  finite-difference  form  as: 
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where  p"  -  pressure  at  the  center  of  the  zone, 


.n+1/2 


velocity  of  the  i^^  grid  point, 
displacement  of  the  i^^  grid  point 


and  dt 


time  step 


The  velocity  and  displacement  of  the  free  boundary  before  averaging  are 


, n+1/2 


(3-5) 


and  d" 


+  u7^1/2 


, respectively. 


(3-6) 


The  following  diagrams  show  the  cycle-by-cycle  cancellation  process 
of  the  ISB  method  as  the  pressure  pulse  propagates  through  the  grid.  In 
these  diagrams,  0  and  u  denote  displacement  and  velocity  of  the  grid 
point,  respectively,  p  denotes  pressure  at  the  center  of  the  zone,  and  U 
denotes  the  velocity  behind  the  step  pulse,  (Pi-Po)/(/>c) .  The  values  of 
0,  u,  and  p  for  each  grid  or  zone  are  computed  according  to  Equations  (3- 
2)  through  (3-6).  The  maximum  stable  time  step  (dt  -  L/c,  c=sound  speed) 
is  assumed  for  each  cycle. 


The  cycle-by-cycle  response  of  the  grids  is  as  follows: 


•  At  cycle  0,  the  grids  are  quiescent,  except  for  the  step  change 
in  pressure  and  velocity  on  the  left  boundary. 


At  cycle  1,  the  pulse  has  propagated  one  zone  width  into  the 
grids.  The  response  in  both  grids  is  identical  because  the 
pulses  have  not  arrived  at  the  right  boundaries.  Averaging  is 
not  shown  because  the  states  are  unchanged. 


At  cycle  2,  the  pulse  has  propagated  two  zone  widths  into  the 
grids.  The  response  is  still  identical  and  averaging  is  not 
shown . 
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At  cycle  3,  the  pulses  first  reflect  off  the  fixed  and  free 
boundaries.  Averaging  produces  perfect  cancellation,  in  the 
sense  that  only  the  incident  wave  remains  after  averaging. 


At  cycle  4,  after  the  second  reflection  from  the  fixed  and  free 
boundaries,  the  averaging  produces  pressures  that  are  different 
than  the  pressure,  pj,  in  the  incident  pulse. 


Initially  at  cycle  0: 
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At  cycle  4  before  averaaina 
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At  cycle  4  after  averaaina 


D-4Udt 
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The  above  analysis  shows  perfect  cancellation  when  the  wave  first 
arrives  at  the  fixed  and  free  boundaries  at  cycle  3.  However,  at  the  next 
cycle,  the  cancellation  error  is  substantial;  the  magnitude  of  the  reflec¬ 
ted  wave  is  50%  of  the  amplitude  of  the  incident  wave.  This  error  will 
persist  if  more  zones  are  used  in  this  analysis.  The  reflected  wave  will 
oscillate  about  zero  stress,  and  the  magnitude  of  the  reflected  wave  can 
be  up  to  50%  of  the  magnitude  of  the  incident  wave.  This  error  will  also 
occur  with  more  complex  pulse  shapes  because  complex  pulses  can  be  genera¬ 
ted  from  superposition  of  simple  step  pulses  (for  small  strain  problems 
with  linear  elastic  materials). 

An  analysis  of  the  simple  analytic  model  shows  that  the  cancellation 
error  occurs  because  the  velocity  of  the  free  boundary  is  not  reset  pro¬ 
perly  after  the  averaging  process. 
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3.2  VELOCITY  CORRECTION  FOR  MAXIMUM  STABLE  TIME  STEP 


The  analytic  model  discussed  in  Section  3.1  shows  that  the  ISB  method 
can  completely  cancel  the  reflected  wave  at  the  fourth  cycle  if  the  velo¬ 
city  at  the  free  boundary  is  reset  to  zero.  The  analytic  model  also 
implies  that  only  one  boundary  zone  is  required  for  cancelling  the 
reflected  waves  although  the  original  ISB  method  averages  three  to  four 
zones  adjacent  to  the  boundaries. 

A  one-dimensional  computational  test  was  performed  with  STEALTH  to 
verify  the  ISB  method  with  the  velocity  at  the  free  boundary  reset  to  zero 
at  the  beginning  of  each  cycle.  The  calculation  is  performed  at  the 
maximum  stable  time  step.  The  length  of  the  computational  grid  is  42  cm, 
excluding  the  one-zone  buffer  region  that  is  1  cm  wide.  The  parameters 
for  the  linear  elastic  material  and  for  the  pressure  boundary  are  identi¬ 
cal  to  the  parameters  for  the  previous  calculations. 

Figure  3-1  presents  the  pressure-time  history  plot  for  a  point  near 
the  transmitting  boundary  (39  cm  from  the  pressure  boundary  or  3  cm  from 
the  inner  boundary  of  the  buffer  zone)  and  the  pressure-distance  plot  for 
the  reflected  wave.  Figure  3-2  presents  the  pressure-distance  plot  of  the 
reflected  wave  in  an  enlarged  scale  at  500  ns,  when  the  reflected  wave  is 
at  the  middle  of  the  main  grid.  These  results  show  that  the  amplitude  of 
the  wave  reflected  from  the  transmitting  boundary  is  0.075  bar,  which  is 
0.75%  of  the  amplitude  of  the  incident  wave.  The  amplitude  of  the 
reflected  wave  is  insignificant  because  the  noise  level  of  a  wave 
reflected  from  a  free  boundary  is  on  the  order  of  0.4%  of  the  incident 
wave. 

3.3  FIRST-ORDER  CORRECTION  FOR  LESS  THAN  MAXIMUM  STABLE  TIME  STEP 

Computational  tests  showed  that  the  cancellation  error  of  the  ISB 
method  is  still  substantial  when  the  time  step  is  less  than  the  maximum 
stable  time  step,  even  though  the  velocity  at  the  free  boundary  is  reset 


to  zero  at  the  beginning  of  each  cycle.  Therefore,  the  velocity  cor¬ 
rection  for  the  free  boundary  has  to  be  modified  for  a  time  step  less  than 
the  maximum  stable  time  step. 

A  first-order  correction  for  the  velocity  at  the  free  boundary  has 
been  developed  for  a  time  step  less  than  the  maximum  stable  time  step. 
Identical  expressions  for  this  first-order  correction  can  be  derived  from 
two  different  approaches.  The  first  approach  is  based  on  a  comparison  of 
the  numerical  solutions  for  a  computational  grid  with  ISB  and  for  a  compu¬ 
tational  grid  with  a  boundary  at  infinity.  The  second  approach  is  based 
on  an  approximation  to  the  stress  at  the  boundary  of  the  computational 
grid. 

First  Approach 

The  first-order  correction  for  the  velocity  at  the  free  boundary  can 
be  derived  by  considering  a  linear  pressure  (or  velocity)  ramp  propagating 
through  a  computational  grid  with  ISB  and  through  a  computational  grid 
with  the  boundary  at  infinity.  The  computational  grid  assumes  equal 
zones.  The  initial  n  cycles  are  computed  at  the  maximum  stable  (and 
constant)  time  step  until  the  leading  edge  of  the  ramp  arrives  at  the  ISB; 
the  subsequent  cycles  are  computed  at  time  steps  less  than  the  maximum 
stable  time  step.  Then  the  velocities  (u)  and  stresses  (a)  at  cycle  n  for 
the  infinite  grid  can  be  expressed  as: 


"  -|r  '‘‘mw 


(n  +  1  -  i)  do. 


(1  ^  i  <  n) 


(1  <  i  <  n) 


(3-7) 


(3-8) 


where  da  ■  stress  increment  per  cycle  at  the  boundary, 
L  -  length  of  a  zone. 


and  dt. 


maximum  stable  time  step. 


At  cycle  n  +  1,  the  velocity  of  the  grid  point  of  the  infinite  grid  is 


(3-9) 


raw 


where  f  =  dt/dt^jax  =  the  time  step,  safety  factor,  and  dt  is  the  time  step 

for  cycle  n  +  1. 

The  velocities  and  stresses  of  the  ISB  grid  are  also  given  by  Equations 
(3-7)  and  (3-8),  except  that  i  is  less  than  or  equal  to  k,  where  k  is  the 
total  number  of  grid  points  in  the  main  grid  and  buffer.  If  the  velocity 
at  the  free  boundary  is  reset  to  at  the  beginning  of  the  cycle  n+1, 

the  velocity  at  the  free  boundary  after  averaging  is  given  by: 


n+1+1/2 


I  [“? 


n+1/2 


6a 

0.5  p  L 


(n  +  1  -  k)  f  dt  ]  (3-10) 


Note  that  the  grid  point  at  the  free  boundary  has  only  half  a  zone  width 
between  the  stress  at  the  center  of  the  zone  and  the  stress  on  the  free 
boundary.  Since  the  velocity  at  the  free  boundary  after  cancellation 
must  be  identical  to  the  velocity  of  the  k^^  grid,  can  be 

evaluated  by  equating  Equations  (2-9)  and  (2-10).  Then  the  first  order 
correction  for  the  velocity  of  the  free  boundary,  ,  can  be 

expressed  as: 


jn+1/2  _  2  ^n.1/2  .  ^ 

=  2(1  -  f) 


(n  -  k)  f  dt 


(3-11) 


If  the  initial  n  cycles  are  not  computed  at  the  maximum  stable  time 
step  and  the  velocity  at  the  free  boundary  is  reset  at  the  beginning  of 
the  cycle  according  to  Equation  (3-11),  the  velocity  at  the  free  boundary 
at  cycle  n  can  be  expressed  as: 
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Second  Approach 


The  second  approach  to  derive  the  first-order  correction  for  the 
velocity  at  the  free  boundary  is  based  on  an  estimate  of  the  stress  in  a 
phantom  zone  at  the  boundary.  This  phantom  zone  has  an  identical  zone 
length  as  the  boundary  zone.  The  new  stress  of  the  phantom  zone  is 
estimated  by  a  linear  interpolation  of  the  stresses  in  the  phantom 
boundary  zone  and  in  the  buffer  zone. 

A  schematic  diagram  for  this  linear  interpolation  is  presented  in 
Figure  3-3.  For  a  linear  elastic  problem,  the  wave  propagates  a  length  of 
cdt  (cdt  -  ft)  per  computational  cycle,  where  c  is  the  sound  speed,  dt  is 
the  computational  time  step,  f  is  the  time  step  safety  factor,  and  L  is 
the  zone  length.  If  constant  time  steps  are  assumed  for  every  cycle,  the 
stress  at  the  phantom  zone,  oU  ,  and  the  velocity  at  the  free  boundary, 


.n+l/2 


,  can  be  expressed  as; 
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Note  that  the  boundary  zone  is  treated  as  an  interior  zone  which  has  a 
fu11  zone  mass  and  width  rather  than  half  width  values.  This  is 
equivalent  to  adding  the  stress  given  by  Equation  (3-13)  to  the  free 
boundary  and  doubling  the  velocity  at  the  free  boundary  at  the  beginning 
of  each  cycle. 

Since  Equation  (3-14)  is  identical  to  Equation  (3-12),  the  velocity 
corrections  at  the  free  boundary  given  by  Equation  (3-11)  are  identical 
for  both  approaches.  Note  that  for  the  maximum  stable  time  step  (f>l)  the 
velocity  correction  given  by  Equation  (3-11)  is  equivalent  to  resetting 


a 


X  -  x,^ 

Note  X  -  - ^ —  , 

a  -  stress, 

X  -  position, 

X|^  -  position  of  the  center  of  the  boundary  zone, 
L  •  one  zone  length,  and 
f  -  time  step  safety  factor. 


Figure  3-3  A  Schematic  Diagram  for  Determining  the  Stress  of  the 
Phantom  Zone  by  a  Linear  Interpolation. 
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boundary  velocity  to  zero.  This  is  consistent  with  the  velocity  correct¬ 
ion  discussed  in  Section  3.2  for  maximum  stable  time  step.  Similarly,  the 
displacements  and  stresses  can  be  shown  to  be  identical  for  both 
approaches. 

A  one-dimensional  test  calculation  was  performed  with  STEALTH.  This 
calculation  is  identical  to  that  discussed  in  Section  3.2  except  that  the 
time  step  was  less  than  the  maximum  stable  time  step  and  the  velocity  at 
the  free  boundary  was  reset  at  the  beginning  of  each  cycle  according  to 
Equation  (3-11). 

Figures  3-4  through  3-6  present  the  pressure-time  history  plots  of 
the  incident  wave  and  the  pressure-distance  plots  of  the  reflected  wave 
for  a  time-step  safety  factor  of  0.95,  0.67,  and  0.50,  respectively. 
Figures  3-7  through  3-9  present,  in  an  enlarged  scale,  the  pressure- 
distance  plots  of  the  reflected  wave  for  the  same  time-step  safety  fact¬ 
ors.  The  pressure-distance  plot  is  taken  at  500  /is,  when  the  wave  front 
of  the  reflected  wave  is  at  the  middle  of  the  main  grid.  These  results 
show  that  the  amplitudes  of  the  reflected  waves  are  0.6,  1.8,  and  2.9%  of 
the  amplitude  of  the  incident  wave  for  a  safety  factor  of  0.95,  0.67,  and 
0.50,  respectively.  Note  that  these  errors  are  much  lower  than  the 
errors  for  the  original  ISB  method.  The  cancellation  errors  for  the 
original  ISB  method  are  on  the  order  of  5  to  12%  of  the  amplitude  of  the 
incident  wave. 

The  residual  cancellation  errors  for  the  ISB  method  with  velocity 
correction  appear  to  be  caused  by  a  second-order  effect,  related  to  the 
rate  of  change  of  pressure.  This  effect  probably  occurs  because  a  time 
step  less  than  the  maximum  stable  time  step  causes  nonlinear  numerical 
diffusion,  while  the  velocity  correction  is  based  on  a  linear  interpola¬ 
tion  for  the  stress  at  the  phantom  boundary  zone. 
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Figure  3-9  Pressure-distance  Plot  at  500  ys  for  th 
for  the  IS6  with  the  First-order  Correc 
Time  Step  Safety  Factor  is  0.50. 


§ 


3.4  SECOND-ORDER  CORRECTION  FOR  LESS  THAN  MAXIMUM  STABLE  TIME  STEP 


3 


Even  though  the  cancellation  errors  of  the  ISB  method  for  a  time  step 
less  than  the  maximum  stable  time  step  can  be  reduced  by  a  first-order 
correction  for  the  velocity  at  the  free  boundary,  the  level  of  the  errors 
is  still  on  the  order  of  0.6%  to  2.9%  of  the  amplitude  of  the  incident 
wave.  A  second-order  correction  for  the  ISB  method  has  been  developed  for 
reducing  the  cancellation  errors  to  the  order  of  the  computational  noise 
level,  which  is  about  0.4%  of  the  amplitude  of  the  incident  wave.  The 
derivation  of  this  second-order  correction  is  based  on  the  extension  of 
the  second  approach  discussed  in  Section  3.3;  the  stress  at  the  phantom 
boundary  zone  is  estimated  by  a  quadratic  fit  rather  than  a  linear  inter¬ 
polation. 

Figure  3-10  presents  a  schematic  diagram  for  estimating  the  phantom 
boundary  zone  by  a  quadratic  fit.  Assuming  the  stress  at  the  phantom 
boundary  zone  satisfies  the  following  quadratic  equation: 

o^  »  a  X  +  bx  +  c  , 

where  =  stress  at  the  boundary  phantom  zone. 


X  =  (X  -  Xj^)  /  L, 

L  =  one  zone  length, 

X  =  position  of  the  center  of  the  phantom  boundary  zone, 

X|^  =  position  of  the  center  of  the  zone  adjacent  to  the  boundary. 

The  coefficients  a,  b,  and  c  are  determined  by  the  following 
conditions: 

at  X  *  0,  » 

at  X  ■  f,  <7^  -  a||'}  and 

at  X  «  2f,  -  a||’? 

where  ajj  is  the  stress  of  the  boundary  zone  at  cycle  n,  and  f  is  the  time 
step  safety  factor. 
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The  stress  in  the  phantom  boundary  zone,  which  is  located  at  x  =  1,  is 
given  by: 


2f2  -  3f  +  1 
2f2 


^  ^  j. 

\  ^ 


2f  -  1  n-1  ^  1 _ 

f2  2,2 


f  n-2 
“  ‘'k 


(3-15) 


Again,  to  calculate  the  velocity  at  the  boundary,  the  boundary  zone  must 
be  treated  as  an  interior  zone  which  has  a  full  zone  width  and  mass.  This 
is  equivalent  to  adding  the  stress  given  by  Equation  (3-15)  to  the  free 
boundary  and  doubling  the  velocity  at  the  free  boundary  at  the  beginning 
of  each  cycle.  Then  the  velocity  at  the  free  boundary  after  cancellation 
can  be  expressed  as: 


n+1/2  _  n-1/2 
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free  free 
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(3-16) 


For  a  time-step  safety  factor  of  one.  Equations  (3-12),  (3-14),  and  (3-16) 
are  identical  and  are  equivalent  to  resetting  the  velocity  of  the  free 
boundary  to  zero  at  the  beginning  of  each  cycle.  This  is  consistent  with 
the  velocity  correction  for  the  maximum  stable  time  step  (Section  3.2). 


One-dimensional  test  calculations  identical  to  those  discussed  in 
Section  3.2  were  performed  with  STEALTH,  except  that  the  time  step  was 
less  than  the  maximum  stable  time  step,  the  velocity  at  the  free  boundary 
was  reset  by  doubling  it  at  the  beginning  of  each  cycle,  and  the  stress 
given  by  Equation  (3-14)  was  added  to  the  free  boundary. 

Figures  3-11  through  3-13  present  the  pressure-time  history  plots  of 
the  incident  wave  and  the  pressure-distance  plots  of  the  reflected  wave 
for  a  time-step  safety  factor  of  0,95,  0.67,  and  0.5,  respectively. 
Figures  3-14  through  3-16  present,  in  an  enlarged  scale,  the  pressure- 
distance  plots  of  the  reflected  wave  for  the  same  time-step  safety  fact¬ 
ors.  These  pressure-distance  plots  are  taken  at  500  /is  when  the  wave 
front  of  the  reflected  wave  is  at  the  middle  of  the  main  grid.  The  resu¬ 
lts  of  the  calculations  show  that  the  amplitude  of  the  reflected  waves  is 
less  than  0.6%  of  the  amplitude  of  the  incident  wave  for  all  safety  fact¬ 
ors.  This  level  of  error  is  insignificant  because  the  computational 
noise  level  of  a  wave  reflected  from  a  free  boundary  is  on  the  order  of 

0.4%  of  the  incident  wave. 
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SECOND-ORDER  CORRECTION. 


(ueqw)  ajnssaud 


Figure  3-11  Pressure  Plots  for  the  IS8  with  the  Sec 
Time  Step  safety  Factor  is  0.95. 
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3-13  Pressure  Plots  for  the  IS8  with  the  Second-order  Correction. 


Pressure  (Mbar) 


STEALTH  ID 

SECOND-ORDER  CORRECTION.  F  -  0.95 


l.OOE-07 


5.00E-08 


0,00E+00 


-5.00E-0a 


-l.OOE-07 


0.00  ib.O  20.0  30.0  40.0 

Distance  (cm) 


Figure  3-14  Pressure-distance  Plot  at  500  us  for  the  Reflected  Wave 
for  the  ISB  with  the  Second-order  Correction. 

Time  Step  Safety  Factor  is  0.95. 
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SECOND-ORDER  CORRECTION.  F  -  0.67 


J  Figure  3-15  Pressure-distance  Plot  at  500  ms  for  the  Reflected  Wave 

for  the  ISB  with  the  Second-order  Correction. 

*  Time  Step  Safety  Factor  Is  0.67. 
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Figure  3-16  Pressure-distance  Plot  at  500  ps  for  the  Reflected  Wave 
for  the  ISB  with  the  Second-order  Correction. 

Time  Step  Safety  Factor  is  0.50. 
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4.0  CONCLUSIONS 


A  transmitting  or  energy  absorbing  boundary  for  one-dimensional 
calculations  of  nonlinear  system  response  in  an  infinite  domain  has  been 
developed.  These  transmitting  boundaries  are  based  on  the  incremental 
superposition  boundary  developed  by  Cundall  et  al .  A  first-  or 
second-order  correction  to  the  ISB  greatly  reduces  the  magnitude  of  the 
reflected  wave  for  all  time  step  safety  factors.  In  addition,  only  one 
zone  in  the  buffer  region  is  required  for  ISB  with  the  first-  or  the 
second-order  correction,  whereas  the  original  ISB  requires  a  buffer  region 
of  three  to  four  zones. 

Calculations  with  a  cosine  pressure  pulse  propagating  through  a 
one-dimensional  grid  were  performed  for  the  original  ISB  and  for  the 
modified  ISB  with  the  first-  and  second-order  corrections.  These  calcula¬ 
tions  were  performed  with  a  Lagrangian,  explicit-in-time,  finHe-differ- 
ence  code  STEALTH.  The  results  of  these  calculations  for  the  magnitude  of 
the  reflected  wave  are  summarized  in  Table  4-1.  These  results  indicate 
that  the  ISB  with  the  first-  or  second-order  corrections  is  more  effective 
in  cancelling  the  reflected  waves  than  the  original  ISB.  In  particular, 
the  second-order  correction  reduces  the  level  of  the  cancellation  error  to 
the  order  of  the  computational  noise,  which  is  about  0.4%  of  the  amplitude 
of  the  incident  wave  for  a  wave  reflected  from  a  free  boundary. 

Although  the  level  of  cancellation  is  good,  the  presence  of  the  time 
step  safety  factor,  f,  in  the  corrections  is  undesirable  for  multi -dimen¬ 
sional  calculations.  In  two-dimensional  calculations,  the  effective 
safety  factor  for  shear  and  dilation  waves  is  different  because  the  sound 
speed  of  both  waves  is  different.  The  presence  of  the  safety  factor  then 
implies  that  a  simple  velocity  correction  factor  will  peroform  equally 
well  for  the  two  different  wave  types.  It  may  be  possible  to  cancel  the 
waves  by  adjusting  the  shear  stress  on  the  boundary.  In  addition  the 
reflected  waves  may  also  be  dispersed  by  the  two-dimensional  reflection 
process.  Further  analysis  is  necessary  to  extend  the  modified  ISB  to  two- 
dimensional  calculations. 


Table  4-1  Comparison  of  the  Magnitude  of  the  Reflected  Wave  for 
the  Original  and  Modified  ISB 


Ratio  of  the  Magnitude  of  the  Reflected 
to  the  Incident  Wave 

Time  Step  Safety  Factors 

0.95 

0.67 

0.50 

The  Original  ISB 

5% 

5% 

12% 

ISB  with  First-order 
Correction 

0.6% 

1.8% 

2.9% 

ISB  with  Second-order 
Correction 


0.6% 


0.5% 


0.5% 
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